Use this template to complete your project throughout the course. Your Final Project presentation will be based on the contents of this document. Replace the title/name above and text below with your own, but keep the headers.

Overview

The topic of this project focuses on bulk RNA-sequencing of lung alveolar type II cells from hyperoxia-injuryed mice. The goal is to see what is the difference in gene expression between alveolar type II cells derived from alveolar type I cells and the rest of the alveolar type II cells. The unpublished sequencing data used in this project comes from Dr. David B. Frank's lab which I currently work in. I first talked with my PI Dr. David B Frank who is a pediatrian at Children's Hospital of Philadelphia focusing on pulmonary diseases to have a better understanding of the motivation behind this preliminary bulk RNA-sequencing and how the results of this experiment can be applied in the medical field. Then I talked with Dr. Prashant Chandrasekaran who is a PhD in Dr. Frank lab with a lot of wet lab experience with mice to learn the details of the experimental procedures. Finally, I talked with Dr. Mengyuan Kan, a postdoctoral in the Department of Biostatistics, Epidemiology and Informatics, at the University of Pennsylvania, working with Dr. Blanca E. Hime, who guides me to analyze the bulk RNA-sequencing data and obtain differential gene expression results. The link to my final project GitHun repository is below:

https://github.com/wenhongbo18/BMIN503_Final_Project

Introduction

Bronchopulmonary dysplasia (BPD) is a chronic lung disease common in premature infants caused by aggressive mechanical ventilator approach on a relatively mature lung lacking surfactant (Davidson et al., 2017; Thebaud et al., 2019). The additional pressure from the ventilation will damage the alveoli, where the gas exchange takes place. Short-term symptoms of BPD include rapid, labored breathing, difficult feeding, repeated lung infections and can eventually lead to asthma-like symptoms, exercise intolerance, and pulmonary arterial hypertension. Almost 50,000 infants were born with less than 28 weeks’ gestation each year in the USA, and approximately 35% of these children develop BPD (Thebaud et al., 2019). Due to the lack of understanding on the mechanisms of BPD, choices of prevention and treatment remain limited and inefficient. Lung alveoli are composed of two types of cells: elongated alveolar type I (AT1) cells which are the site of gas exchange and cuboidal alveolar type II (AT2) cells whose major role is to produce surfactant which lowers the surface tension of alveoli and makes sure they do not collpase when expanding. Previous studies have demonstrated the plasticity of AT2 progenitor cells and how they can differentiate and proliferate to heal the injuried alveoli in adult mice (Barkauskas et al., 2013). Compared to AT2 cells, AT1 cells were previously believed to lack plasticity and cannot change their cell fates. However, a recent study shows that upon acute neonatal lung injury, AT1 cells reprogram into AT2 cells, promoting alveolar regeneration. Considering these background information, our research group asks the question what signaling pathways govern the plasticity of AT1 cells. The answer to this question will have great application in treating BPD. Researchers can find ways to upregulate or downregulate these signaling pathways to increase plasticity of AT1 cells, thereby promoting the healing of alveolar damage caused by BPD.

An interdisciplinary approach is needed to solve this problem. A strong biology and chemistry background is needed to design the experiment, choose the right mice model, and obtain the specific types of cells. Afterwards, biomedical informatics knowledge comes in to help organize the data, perform statistical analysis, and locate the genes that have the highest difference in expression between two population of the cells. Then, a good interpretation of the data requires a comprehensive understanding of genetics. After finding relevant genes, the search of existing medications that target related signaling pathway and the determination of the effectiveness of such medications need perspectives from the clinical side.

Methods

Describe the data used and general methodological approach. Subsequently, incorporate full R code necessary to retrieve and clean data, and perform analysis. Be sure to include a description of code so that others (including your future self) can understand what you are doing and why.

The mouse model used in this experiment was the HopXCre R26Rtdtomato SPC-GFP mice. Due to the gentic design, all the AT2 cells were marked with a green color by the green fluroscent protein (GFP), but only the AT1-derived AT2 cells also expressed a red color (tdtomato) from the Cre-LoxP lineage tracing design. Because ths is only a preliminary experiment conducted in my lab, only two mice were used named 927_1 and 980_3. Since we were interested in exploring the AT1 and AT2 cells fate after injury, the two mice were put in the hyperoxia chamber at 90% of oxygen level for four days to mimic the conditions of BPD. The mice were then transported to live in normal condition until postnatal day 42. The lungs were then harvested, single cell suspension was achieved, and samples were submitted to flow cytometry after appropriate antibody staining.Two cell populations were harvested through flow cytometry: AT1-derived AT2 cells (GFP-tdtomato) and other AT2 cells (GFP). Bulk RNA-sequencing is performed on these two sets of cells to see what difference in gene expression explains this phenotype. The two different groups of cells underwent cell lysis, mRNA extraction, mRNA fragmentation, reverse transcription to cDNA, PCR amplification, ligation of sequence adaptors, and size selection. Afterwards, the samples were put into the flow cells of Illumina sequencer to obtain the sequence of all the cDNAs, and the data was stored in a hard drive.

All the data was processed on a remote server at Dr. Himes lab, so the first step is to transfer all the data into the remote server. A PMACS account was set up for me to use. Linux command is used below.

  1. Login to the server, use the command:

ssh hwen@bhimes.pmacs.upenn.edu

  1. Under the RNASeq folder, create a folder called Lung_AT2:

mkdir -p /projects/RNASeq/Lung_AT2

  1. Transfer all the sequencing data into this folder:

scp -r (full path of the data folder on my computer) hwen@bhimes.pmacs.upenn.edu:/projects/RNASeq/Lung_AT2 df -h

  1. Initialize the python 2.7:

conda activate py2

  1. Create a folder named files and inside that folder create a txt file called Lung_AT2_Phenotype_withoutQC.txt.

The Lung_AT2_Phenotype_withoutQC.txt file has 8 columes:

Sample: Sample name Donor: Names of the mice that the samples were from Tissue: Type of the tissue the cells came from. Here it is alveolar type 2 cells from lung tissues. Treatment: Hyperoxia injury Disease: All mice were healthy beforehand Status: This used for comparison. We have two groups: AT2 cells derived from AT1 cells (AT2_AT1) and other AT2 cells (AT2) R1: directory of the R1 read (5' direction of the top strand) R2: directory of the R2 read (5' direction of the bottom strand)

mkdir -p /projects/RNASeq/Lung_AT2/files

cd /projects/RNASeq/Lung_AT2/files

cat > Lung_AT2_Phenotype_withoutQC.txt

library(readr)
withoutQC <- read_csv("~/Downloads/RNASeq.csv")
head(withoutQC)
## # A tibble: 4 × 8
##   Sample           Donor  Tissue Treatment Disease Status  R1                                                               R2                                                              
##   <chr>            <chr>  <chr>  <chr>     <chr>   <chr>   <chr>                                                            <chr>                                                           
## 1 DF_927_1_GFP_tdt s927_1 AT2    Hyperoxia healthy AT2_AT1 /projects/RNASeq/Lung_AT2/00_fastq/927-1-GFP-TdT_R1_001.fastq.gz /projects/RNASeq/Lung_AT2/00_fastq/927-1-GFP-TdT_R2_001.fastq.gz
## 2 DF_927_1_GFP     s927_1 AT2    Hyperoxia healthy AT2     /projects/RNASeq/Lung_AT2/00_fastq/927-1-GFP_R1_001.fastq.gz     /projects/RNASeq/Lung_AT2/00_fastq/927-1-GFP_R2_001.fastq.gz    
## 3 DF_980_3_GFP_tdt s980_3 AT2    Hyperoxia healthy AT2_AT1 /projects/RNASeq/Lung_AT2/00_fastq/980-3-GFP-TdT_R1_001.fastq.gz /projects/RNASeq/Lung_AT2/00_fastq/980-3-GFP-TdT_R2_001.fastq.gz
## 4 DF_980_3_GFP     s980_3 AT2    Hyperoxia healthy AT2     /projects/RNASeq/Lung_AT2/00_fastq/980-3-GFP_R1_001.fastq.gz     /projects/RNASeq/Lung_AT2/00_fastq/980-3-GFP_R2_001.fastq.gz
  1. Run pipeline scripts/rnaseq align and qc.py to: 1) trim adapter and primer sequences if index information is available, 2) run FastQC for (un)trimmed .fastq files, 3) align reads and quantify reads mapped to genes/transcripts, and 4) obtain various QC metrics from .bam files. At the end of these steps, the readcount file was generated:

mkdir -p /projects/RNASeq/Lung_AT2/scripts/align

cd /projects/RNASeq/Lung_AT2/scripts/align

rnaseq_align_and_qc.py --project_name Lung_AT2 --samples_in /projects/RNASeq/Lung_AT2/files/Lung_AT2_Phenotype_withoutQC.txt --aligner star --ref_genome mm10 --library_type PE --strand nonstrand --path_start /projects/RNASeq/Lung_AT2 --template_dir /projects/pipelines/RNASeq

mk_lsf2bashsub.py --path /projects/RNASeq/Lung_AT2/scripts/align --lsf_dir /projects/RNASeq/Lung_AT2/scripts/align --job_name align --nthread 1 & echo $! > /projects/RNASeq/Lung_AT2/scripts/align/align.jobid

jobid=cat /projects/RNASeq/Lung_AT2/scripts/align/align.jobid ps -A | grep $jobid

  1. create an HTML report of QC and alignment summary statistics for RNA-seq samples:

mkdir -p /projects/RNASeq/Lung_AT2/scripts/qc_report

cd /projects/RNASeq/Lung_AT2/scripts/qc_report

rnaseq_align_and_qc_report.py --project_name Lung_AT2 --samples_in /projects/RNASeq/Lung_AT2/files/Lung_AT2_Phenotype_withoutQC.txt --aligner star --ref_genome mm10 --library_type PE --strand nonstrand --path_start /projects/RNASeq/Lung_AT2 --template_dir /projects/pipelines/RNASeq

bash /projects/RNASeq/Lung_AT2/scripts/qc_report/Lung_AT2_qc.lsf > Lung_AT2_qc.log 2>&1 & echo $! > /projects/RNASeq/Lung_AT2/scripts/qc_report/qc_report.jobid

jobid=cat /projects/RNASeq/Lung_AT2/scripts/qc_report/qc_report.jobid ps -A | grep $jobid

  1. Inside the file folder create two txt files called Lung_AT2_Phenotype_withQC.txt, and Lung_AT2_comp_file.txt.

The Lung_AT2_Phenotype_withoutQC.txt file has 9 columes:

Sample: Sample name Donor: Names of the mice that the samples were from Tissue: Type of the tissue the cells came from. Here it is alveolar type 2 cells from lung tissues. Treatment: Hyperoxia injury Disease: All mice were healthy beforehand Status: This used for comparison. We have two groups: AT2 cells derived from AT1 cells (AT2_AT1) and other AT2 cells (AT2) R1: directory of the R1 read (5' direction of the top strand) R2: directory of the R2 read (5' direction of the bottom strand) QC_Pass: 1 if the sample QC is good and can be included in the following DE analysis; 0 if the sample QC is bad and should be excluded in the following DE analysis. (Note: this is done after the QC report is generated. All the samples have good QCs and can be used in the following DE analysis)

The Lung_AT2_comp_file.txt file has 3 columes:

Condition1: AT2 cells derived from AT1 cells (AT2_AT1) Condition0: other AT2 cells (AT2) Design: paired:Donor. This specifies whether this sample is paired with another sample for comparison. All of the 4 samples used were paired.

mkdir -p /projects/RNASeq/Lung_AT2/files

cd /projects/RNASeq/Lung_AT2/files

cat > Lung_AT2_Phenotype_withQC.txt

cat > Lung_AT2_comp_file.txt

library(readr)
withQC <- read_csv("~/Downloads/withQC.csv")
head(withQC)
## # A tibble: 4 × 9
##   Sample           Donor  Tissue Treatment Disease Status  R1                                                               R2                                                               QC_Pass
##   <chr>            <chr>  <chr>  <chr>     <chr>   <chr>   <chr>                                                            <chr>                                                              <dbl>
## 1 DF_927_1_GFP_tdt s927_1 AT2    Hyperoxia healthy AT2_AT1 /projects/RNASeq/Lung_AT2/00_fastq/927-1-GFP-TdT_R1_001.fastq.gz /projects/RNASeq/Lung_AT2/00_fastq/927-1-GFP-TdT_R2_001.fastq.gz       1
## 2 DF_927_1_GFP     s927_1 AT2    Hyperoxia healthy AT2     /projects/RNASeq/Lung_AT2/00_fastq/927-1-GFP_R1_001.fastq.gz     /projects/RNASeq/Lung_AT2/00_fastq/927-1-GFP_R2_001.fastq.gz           1
## 3 DF_980_3_GFP_tdt s980_3 AT2    Hyperoxia healthy AT2_AT1 /projects/RNASeq/Lung_AT2/00_fastq/980-3-GFP-TdT_R1_001.fastq.gz /projects/RNASeq/Lung_AT2/00_fastq/980-3-GFP-TdT_R2_001.fastq.gz       1
## 4 DF_980_3_GFP     s980_3 AT2    Hyperoxia healthy AT2     /projects/RNASeq/Lung_AT2/00_fastq/980-3-GFP_R1_001.fastq.gz     /projects/RNASeq/Lung_AT2/00_fastq/980-3-GFP_R2_001.fastq.gz           1
comp <- read_csv("~/Downloads/comp.csv")
head(comp)
## # A tibble: 1 × 3
##   Condition1 Condition0 Design      
##   <chr>      <chr>      <chr>       
## 1 AT2_AT1    AT2        paired:Donor
  1. Perform DE analysis and create an HTML report of differential expression summary statistics. At the end of these steps, the readcount file was generated:

mkdir -p /projects/RNASeq/Lung_AT2/scripts/deseq2

cd /projects/RNASeq/Lung_AT2/scripts/deseq2

rnaseq_de_report.py --project_name Lung_AT2 --samples_in /projects/RNASeq/Lung_AT2/files/Lung_AT2_Phenotype_withQC.txt --comp /projects/RNASeq/Lung_AT2/files/Lung_AT2_comp_file.txt --de_package deseq2 --ref_genome mm10 --path_start /projects/RNASeq/Lung_AT2 --template_dir /projects/pipelines/RNASeq

bash /projects/RNASeq/Lung_AT2/scripts/deseq2/Lung_AT2_deseq2.lsf > Lung_AT2_deseq2.log 2>&1 & echo $! > /projects/RNASeq/Lung_AT2/scripts/deseq2/Lung_AT2_deseq2.jobid

jobid=cat /projects/RNASeq/Lung_AT2/scripts/deseq2/Lung_AT2_deseq2.jobid ps -A | grep $jobid

  1. Copy all the files needed from the remote server to my laptop:

scp -r hwen@bhimes.pmacs.upenn.edu:/projects/RNASeq/Lung_AT2/Lung_AT2_deseq2_out/ /Users/tony/Downloads

scp -r hwen@bhimes.pmacs.upenn.edu:/projects/RNASeq/Lung_AT2/scripts /Users/tony/Downloads

scp -r hwen@bhimes.pmacs.upenn.edu:/projects/RNASeq/Lung_AT2/files /Users/tony/Downloads

scp -r hwen@bhimes.pmacs.upenn.edu:/projects/RNASeq/Lung_AT2/Lung_AT2_Alignment_QC_Report_star /Users/tony/Downloads

Results

Quality Control

Project: Lung_AT2

Aligner: STAR (2.7.6a)

Genome: For mouse, the UCSC mm10 assembly available in iGenomes was used.

Informatics tools used:

Sequencing parameters:

For each sample, the following programs were run to generate the data necessary to create this report. Written as for unstranded paired-end data. For single-end reads, R2s and insert size metrics would be omitted.

java -Xmx1024m TrimmomaticPE -phred33 [raw_sample_R1] [raw_sample_R2] [sample_R1] [sample_R1_unpaired] [sample_R2] [sample_R2_unpaired] HEADCROP:[bases to trim, if any] ILLUMINACLIP:[sample_primer_fasta]:2:30:10 MINLEN:50

fastqc [sample_R1] [sample_R2]

cat [sample_R1/R2] | awk '((NR-2)%4==0){read=$1;total++;count[read]++}END{for(read in count){if(count[read]==1){unique++}};print total,unique,unique*100/total}'

The following STAR options were used:

STAR --genomeDir [ref_genome_index] --runThreadN 12 --outReadsUnmapped Fastx --outMultimapperOrder Random --outSAMmultNmax 1 --outFilterIntronMotifs RemoveNoncanonical --outSAMstrandField intronMotif --outSAMtype BAM SortedByCoordinate --readFilesIn [sample_R1] [sample_R2]

Using aligned output files accepted_hits.bam and unmapped.bam:

samtools sort accepted_hits.bam accepted_hits.sorted

samtools index accepted_hits.sorted.bam

samtools idxstats accepted_hits.sorted.bam > accepted_hits.sorted.stats

bamtools stats -in accepted_hits.sorted.bam > accepted_hits.sorted.bamstats

bamtools filter -in accepted_hits.sorted.bam -script cigarN.script | bamtools count

samtools view -c unmapped.bam

java -Xmx2g -jar CollectRnaSeqMetrics.jar REF_FLAT=[ref_flat file] STRAND_SPECIFICITY=NONE INPUT=accepted_hits.bam OUTPUT=RNASeqMetrics

java -Xmx2g -jar CollectInsertSizeMetrics.jar HISTOGRAM_FILE=InsertSizeHist.pdf INPUT=accepted_hits.sorted.bam OUTPUT=InsertSizeMetrics (for paired-end library)

The following R codes generate the readcounts.

## Warning: package 'BiocGenerics' was built under R version 4.0.5
## Warning: package 'GenomeInfoDb' was built under R version 4.0.5

Summary Read Numbers

The number of raw reads correspond to those that passed Casava QC filters, were trimmed to remove adaptors by Trimmomatic, and were aligned by STAR to ref_genome+ERCC transcripts as reported in .info files. Unique read counts were obtained by using awk on trimmed fastq files. FastQC estimates of percentage of sequences remaining after deduplication were retrieved from fastqc_data.txt files. Bamtools statistics were based on sorted and indexed bam files. The mapped reads were those that mapped to reference and were output by STAR to accepted_hits.bam. The unmapped reads were output by STAR to unmapped.bam. Some reads may be mapped to multiple locations in the genome so that the number of total reads reported by bamstats may be greater than the number of raw reads. The Junction spanning reads are computed based on accepted_hits.bam CIGAR entries containing "N." Related text files that were saved:

Lung_AT2 _read_counts.txt

Lung_AT2 _duplicates.txt

Lung_AT2 _unique_counts.txt

Lung_AT2 _bamstats_counts.txt

Total Number of Raw Reads Summary

Read counts are shown by per million reads.

Plot: Percentage of Unique Reads in Original Fastq File

Plot: Sequence Duplication Level

Bamtools Reads Summary

Bamtools Reads Summary As Percentage of Mapped Reads

Percentage of Mapped/Unmapped Reads

Plot: Percentage of Mapped/Unmapped Reads

Plot: Percentage of Junction Spanning Reads Among Mapped Reads

RnaSeqMetrics Summary

The Picard Tools RnaSeqMetrics function computes the number of bases assigned to various classes of RNA. It also computes the coverage of bases across all transcripts (normalized to a same-sized reference). Computations are based on comparison to a refFlat file. Related text files that were saved:

Lung_AT2 _rnaseqmetrics_summary.txt

Lung_AT2 _rnaseqmetrics_hist.txt

The Picard Tools RnaSeqMetrics function computes the number of bases assigned to various classes of RNA. It also computes the coverage of bases across all transcripts (normalized to a same-sized reference). Computations are based on comparison to a refFlat file. Related text files that were saved:

Lung_AT2 _rnaseqmetrics_summary.txt

Lung_AT2 _rnaseqmetrics_hist.txt

Reference Genome Mapped Reads Summary

Plot: Percentages of Total Mapped Bases Mapping to mRNA, Intronic and Intergenic Regions

Plot: Normalized Coverage

InsertSizeMetrics Summary

For paired-end data, the Picard Tools CollectInsertSizeMetrics function was used to compute the distribution of insert sizes in the accepted_hits.bam file and create a histogram. Related text files that were saved:

Lung_AT2 _insertmetrics_summary.txt

Plot: Median of Insert Size

Plot: Insert Size Distribution

Reads per Chromosome

Samtools produces a summary document that includes the number of reads mapped to each chromosome. Related text files that were saved:

Lung_AT2 _counts.txt

Mapped Reads to Reference Genome

Percent of Total Reads Mapped to Reference Genome

Differential Expression Analysis

Reads were aligned to the mm10 assembly using STAR (2.7.6a). The following alignment QC report was produced:

Lung_AT2_QC_RnaSeqReport.html

HTSeq (0.11.3) function htseq-count was used to count reads. Counts for all samples were concatenated into the following text file:

Lung_AT2_htseq_gene.txt

DESeq2 (1.28.1) was used for differential gene expression analaysis, based on the HTSeq counts matrix and the phenotype file provided. Normalized counts from DESeq2 are saved in the following text file:

Lung_AT2_counts_normalized_by_DESeq2.txt

Normalized counts are obtained from DESeq2 function estimateSizeFactors(), which divides counts by the geometric mean across samples; this function does not correct for read length. The normalization method is described in detail here: https://genomebiology.biomedcentral.com/articles/10.1186/gb-2010-11-10-r106

Differential gene expression analysis was done for all comparisons provided in the comparisons file. The following design was used:

design = ~ + Donor + Status

If desired, the design can be modified to include more independent variables. In addition to the partial results displayed in this report, the full set of DESeq2 results for each comparison was saved down in separate text files, with names of the form:

Lung_AT2_CASE_vs_CONTROL_DESeq2_results.txt

where CASE and CONTROL are pairs of conditions specified in the comparisons file.

## [1] "/Users/tony/BMIN503_Final_Project"
## Warning in read.table("/Users/tony/BMIN503_Final_Project/files/Lung_AT2_Phenotype_withQC.txt", : incomplete final line found by readTableHeader on '/Users/tony/BMIN503_Final_Project/files/Lung_AT2_Phenotype_withQC.txt'

Convert KEGG and REACTOME pathway files in .gmt format provided by MSigDB into a pathway list

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in design formula are characters, converting to factors

AT2_AT1 vs. AT2

Samples in this comparison

DE analysis

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in design formula are characters, converting to factors
Comparison Summary
Status Count
AT2 2
AT2_AT1 2

Top 50 genes by p-value

Description of DESEq2 output

Volcano plots

Volcano plot (probes with a q-value <0.05 are present in red)

MA plot

Distribution of adjusted p-values

Dendrogram based on sample distance of regularized log transformed data

Heatmaps for top 30 significant genes

Genes were ranked by adjusted p-values.

Principal component analysis (PCA) Plot based on regularized log transformed data

Compute PCs and variance explained by the first 10 PCs

Variance explained
PC Proportion of Variance (%) Cumulative Proportion of Variance (%)
PC1 70.82 70.82
PC2 19.06 89.88
PC3 10.12 100
PC4 1.077e-28 100

PCA plots are generated using the first two principle components colored by known factors (e.g. Status, Tissue, or Donor)

Dispersion plot

Plot of the maximum Cook's distance per gene over the rank of the Wald statistics for the condition

Boxplots for top 20 differentially expressed genes

Genes were ranked by pvalue. Counts have been normalized by sequencing depth, with pseudocount of 0.5 added to allow for log scale plotting, using DESeq2 function plotCounts().

Favorite genes

Boxplots for user-defined favorite genes if they exist, and show DE results

## None of user-defined genes are found in current DE results.

Gene-set enrichment analysis

View top pathways in fgsea results. Only include the first five leading genes

## 52 main pathways are significant.

Generate barplot if pathways with absoluate NES>=2 or top 10 pathways if no pathways pass the threshold

View leading edges in top pathways. Select top five pathways with positive and negative NES respectively

Housekeeping genes

Counts have been normalized by estimated size factors using DESeq2. Obtain the count matrix using function DESeq2::counts.

The table shows p-values of house-keeping genes for each comparison. Generally, house-keeping gene expressions do not change significantly in different conditions.

Favorite gene expressions in all conditions

R version 4.0.3 (2020-10-10)

**Platform:** x86_64-apple-darwin17.0 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: parallel, stats4, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: knitr(v.1.36), pander(v.0.6.4), fgsea(v.1.16.0), biomaRt(v.2.46.3), DESeq2(v.1.30.1), SummarizedExperiment(v.1.20.0), Biobase(v.2.50.0), MatrixGenerics(v.1.2.1), matrixStats(v.0.61.0), GenomicRanges(v.1.42.0), GenomeInfoDb(v.1.26.7), IRanges(v.2.24.1), S4Vectors(v.0.28.1), BiocGenerics(v.0.36.1), viridis(v.0.6.2), viridisLite(v.0.4.0), gplots(v.3.1.1), ggplot2(v.3.3.5), DT(v.0.20), RColorBrewer(v.1.1-2), tidyr(v.1.1.4) and readr(v.2.1.0)

loaded via a namespace (and not attached): bitops(v.1.0-7), bit64(v.4.0.5), progress(v.1.2.2), httr(v.1.4.2), tools(v.4.0.3), utf8(v.1.2.2), R6(v.2.5.1), KernSmooth(v.2.23-20), DBI(v.1.1.1), colorspace(v.2.0-2), withr(v.2.4.2), prettyunits(v.1.1.1), tidyselect(v.1.1.1), gridExtra(v.2.3), curl(v.4.3.2), bit(v.4.0.4), compiler(v.4.0.3), cli(v.3.1.0), xml2(v.1.3.2), DelayedArray(v.0.16.3), labeling(v.0.4.2), caTools(v.1.18.2), scales(v.1.1.1), genefilter(v.1.72.1), askpass(v.1.1), rappdirs(v.0.3.3), stringr(v.1.4.0), digest(v.0.6.28), rmarkdown(v.2.11), XVector(v.0.30.0), pkgconfig(v.2.0.3), htmltools(v.0.5.2), highr(v.0.9), dbplyr(v.2.1.1), fastmap(v.1.1.0), htmlwidgets(v.1.5.4), rlang(v.0.4.12), rstudioapi(v.0.13), RSQLite(v.2.2.8), farver(v.2.1.0), jquerylib(v.0.1.4), generics(v.0.1.1), jsonlite(v.1.7.2), crosstalk(v.1.2.0), BiocParallel(v.1.24.1), gtools(v.3.9.2), vroom(v.1.5.6), dplyr(v.1.0.7), RCurl(v.1.98-1.5), magrittr(v.2.0.1), GenomeInfoDbData(v.1.2.4), Matrix(v.1.3-4), Rcpp(v.1.0.7), munsell(v.0.5.0), fansi(v.0.5.0), lifecycle(v.1.0.1), stringi(v.1.7.5), yaml(v.2.2.1), zlibbioc(v.1.36.0), BiocFileCache(v.1.14.0), grid(v.4.0.3), blob(v.1.2.2), crayon(v.1.4.2), lattice(v.0.20-45), splines(v.4.0.3), annotate(v.1.68.0), hms(v.1.1.1), locfit(v.1.5-9.4), pillar(v.1.6.4), geneplotter(v.1.68.0), fastmatch(v.1.1-3), XML(v.3.99-0.8), glue(v.1.5.0), evaluate(v.0.14), data.table(v.1.14.2), vctrs(v.0.3.8), tzdb(v.0.2.0), openssl(v.1.4.5), gtable(v.0.3.0), purrr(v.0.3.4), assertthat(v.0.2.1), cachem(v.1.0.6), xfun(v.0.28), xtable(v.1.8-4), survival(v.3.2-13), tibble(v.3.1.6), AnnotationDbi(v.1.52.0), memoise(v.2.0.0) and ellipsis(v.0.3.2)